#!/usr/bin/perl
# 2012/10/17 Output location of target sequence

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [fasta file] -t [recognition sequence]\n";
}

my %opt;
getopts("i:t:", \%opt);
my $FastaFile = $opt{i};
my $recognitionSeq = $opt{t};
$recognitionSeq = uc($recognitionSeq);

print "#number\tchr\tposition\tlength_before\tlength_after\n";
my $fh = IO::File->new($FastaFile) or die "cannot open $FastaFile: $!";
my $id = '';
my $seq = '';
while($_ = $fh->getline()){
	s/\r?\n//;
	if(m/^>(\S+)/){
		$id = $1;
	}else{
		$seq = $_;
		$seq = uc($seq);
		&parseSeq($id, $seq);
	}
}
$fh->close();

sub parseSeq{
	my ($chr, $seq) = @_;
	my $frag_number = 0;
	if($chr eq 'I'){
		$frag_number += 100000;
	}elsif($chr eq 'II'){
		$frag_number += 200000;
	}elsif($chr eq 'III'){
		$frag_number += 300000;
	}elsif($chr eq 'MT'){
		$frag_number += 400000;
	}elsif($chr eq 'MTR'){
		$frag_number += 500000;
	}elsif($chr eq 'AB325691'){
		$frag_number += 600000;
	}else{
		$frag_number += 700000;
	}

	my $pos = 0;
	my $previous = 0;
	my @locationList = (0);
	while($pos != -1){
		$pos = index $seq, $recognitionSeq, $previous;
		if($pos != -1){
			push @locationList, $pos;
		}
		$previous = $pos + 1;
	}
	my $totalLength = length $seq;
	push @locationList, $totalLength;

	for(my $i = 1; $i < @locationList - 1; $i++){
		$frag_number++;
		my $length_before = $locationList[$i] - $locationList[$i-1];
		my $length_after = $locationList[$i+1] - $locationList[$i];
		print "$frag_number\t$chr\t$locationList[$i]\t$length_before\t$length_after\n";
	}
}

